##############################################################################
# This is the simulation script. See the README for more information.
##############################################################################
from __main__ import *

import copy
import operator
import numpy as np
import itertools as it
import time
import mpmath
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib import font_manager
from skellam import skcdf
from skellam import skpmf


###############################################################################
#Global variables
###############################################################################
CANDS = ["gross", "claus", "begich", "palin", "peltola"]

SEED = 20220909 #The date that wthis script was first drafted

VERBOSE = False
ERROR_CHECKS = False

NUM_RUNS = 1 #Number of runs for the SMD vs IRV pivot probability simulations

NO_SWITCHES = {"base": ["foo"], "alt": ["bar"]} #Switches that won't print

#Uncomment the font family setting to have the same font as I use in the paper.
# It is commented out in this version because CodeOcean does not support it
#plt.rcParams['font.family'] = 'Ubuntu Mono'
plt.rcParams['font.size'] = '22'
dirPivotProbs = {}
indirPivotProbs = {}


###############################################################################
# Functions
###############################################################################
def RepList(theList, n):
    reppedList = [_ for _ in it.repeat(theList, n)]
    return(reppedList)

#Function to count votes given a list of ballots and of dropped candidates
def GetVotes(ballots, cands, dropped):
    currCands = [c for c in cands if c not in dropped]
    votes = {c:0 for c in currCands}
    currBallots = copy.deepcopy(ballots)
    for ballot in currBallots:
        topCand = ballot[0]
        while topCand in dropped:
            ballot = [b for b in ballot if b != topCand]
            topCand = ballot[0]
        votes[topCand] += 1
    return(votes)

def IRV(cands, ballots):
    _ballots = copy.deepcopy(ballots)
    _irvRes = {c:0 for c in cands}
    for ballot in _ballots:
        _irvRes[ballot[0]] += 1
    print(str(_irvRes))
    #Iteratively reallocate votes until we have only one winner
    while len(_irvRes.keys()) > 1:
        #Identify the lowest-performing candidate, with order-based
        # tie-breaking, and drop it from the dictionary of results
        _lowestCand = min(_irvRes.items(),key=operator.itemgetter(1))[0]
        _irvRes.pop(_lowestCand)
        #Then loop over ballots, and for each ballot on which the dropped 
        # candidate was at the top, identify to whom that ballot transfers
        for i in range(len(_ballots)):
            if _ballots[i][0] == _lowestCand:
                _irvRes[_ballots[i][1]] += 1
            if _lowestCand in _ballots[i]:
                _ballots[i] = [c for c in _ballots[i] if c != _lowestCand]
        print(str(_irvRes))
    return(_irvRes)

def CalcDropProbs(sequences, ballots, cands):
    #Globally calculate the probabilities of candidates dropping in every order
    # to then use in utility computations. What we want is a dictionary where
    # we can input any order of candidates dropping of any length, and the
    # dictionary returns the probability of that drop order occurring
    _all_drop_probs = {}
    for order in sequences:
        dropped = []
        orderProb = 1
        orderLabel = ListToStr(order)
        for dropCand in order:
            votes = GetVotes(HELD_BALLOTS, CANDS, dropped)
            #We need to compare this candidate to every candidate that has not
            #yet been dropped
            remCands = [_ for _ in CANDS if _ not in dropped and _ != dropCand]
            for remCand in remCands:
                orderProb *= (skcdf(-1, votes[dropCand], votes[remCand]) +\
                              .5*skpmf(0, votes[dropCand], votes[remCand]))
                #print(order, dropCand, remCand, orderProb)
            dropped.append(dropCand)
        _all_drop_probs[orderLabel] = orderProb
    return(_all_drop_probs)

#String and destring lists for using them as dictionary keys
def ListToStr(theList):
    listStr = ''
    for i in theList:
        listStr += str(i)
    return(listStr)

def StrToList(theString, toInt):
    strList = []
    for i in theString:
        if toInt:
            i = int(i)
        strList.append(i)
    return(strList)

def CalcDirectPiv(cands, ballot, all_ballots, all_drop_probs, probs_to_print):
    _dirPivotProbs = {}

    #To get direct pivotality for candidates in ballot positions beyond the
    # first position, condition on the candidate that you have put in your 
    # first ballot position being dropped in each possible round, and then
    # check the pivotal probability of the candidate in the next position
    # from there.
    ballotLabel = ListToStr(ballot)
    _dirPivotProbs[ballotLabel] = 0

    #Direct pivotal probabilities
    for i in range(len(ballot)):
        c = ballot[i]

        #If this is not the first candidate on the ballot, we need all the
        # other candidates to have been dropped. So we only care about drop
        # sequences which have the required candidates among them (in the
        # case of they first candidate on the ballot, that list is empty)
        prevDrop = [d for d in ballot[:i]]
        
        #We are interested in all drop sequences in which c is a final
        # contender, that is, all length M-1 drop sequences, which:
        # 1) do not include c
        # 2) must include every candidate above c on the ballot
        # 3) have c in the last position
        dropped = [d for d in cands if d != c]
        dropSeqs = [list(d) for d in it.permutations(dropped)]

        for seq in dropSeqs:
            #At issue is the probability that the final contest is tied.
            # (note: currently doing redundant computations because drop
            # order doesn't matter)
            prevDropSeq = seq[:-1]

            #Each drop sequence has a pre-computed probability of happening,
            # so we need the probability of a further drop sequence 
            # conditioned on this sequence of events
            dropSeqProb = all_drop_probs[ListToStr(prevDropSeq)]
            
            #The candidates which still have votes are c and its direct
            # opponent
            votes = GetVotes(all_ballots, cands, prevDropSeq)

            #Now that we have the probability of this sequence of drops
            # occurring, we need only to know the probability that a vote
            # for the candidate under consideration will be pivotal at
            # the end of this chain
            k = seq[-1]
            if votes[c] > 0 and votes[k] > 0:
                breakProb = dropSeqProb * skpmf(0, votes[c], votes[k])
                makeProb = dropSeqProb * skpmf(-1, votes[c], votes[k])
            else:
                (breakProb, makeProb) = (0, 0)
            fullPivotProb = breakProb + 0.5*makeProb
            _dirPivotProbs[ballotLabel] += fullPivotProb
            if seq in probs_to_print:
                print(f"{prevDropSeq} have been dropped with prob "+\
                      f"{dropSeqProb}, {c} has"+\
                      f" {votes[c]} votes, {k} has {votes[k]} votes."+\
                      f" Tie break prob of {c} vs. {k} is {breakProb} and"+\
                      f" tie make prob is {makeProb}, with full pivot prob"+\
                      f" {fullPivotProb}")
    return(_dirPivotProbs)
        

def CalcIndirectPiv(cands,ballot,all_ballots,all_drop_probs,switch_to_print):
    indirPivotProbs = {}
    '''
    At its heart indirect pivotality is the probability of bouncing from a
    chain with one winner to a chain with a different winner. Without any
    simplifying theorems, this means that the procedure is to take any
    pair of drop sequences and consider the chances that a vote for a
    candidate who does not win in those sequences causes the result to change
    from one of those chains to the other one
    '''
    ballotLabel = ListToStr(ballot)
    indirPivotProbs[ballotLabel] = 0
    
    #List all of the sequences that could occur, before boiling them down to
    # only the relevant ones
    allSeqs = [list(p) for p in it.permutations(cands)]
    
    #We only consider candidates who are early enough on the ballot that at
    # least 4 candidates follow them
    for i in range(len(ballot[:-4])):
        c = ballot[i]
        
        otherCands = [j for j in cands if j != c]
        
        #By a theorem, we can only be indirectly pivotal in rounds where there
        # are more than 4 candidates remaining
        allMainSeqs = [s for s in allSeqs if c not in s[-4:]]
        reqDropCands = ballot[:i]
        #Candidates that appear before c on the ballot must have been dropped
        # before c
        mainSeqs = []
        for seq in allMainSeqs:
            keepThisSeq = True
            for req in reqDropCands:
                if seq.index(req) > seq.index(c):
                    keepThisSeq = False
            if keepThisSeq:
                mainSeqs.append(seq)
        
        #What is required for one sequence to morph into another?
        # 1. Every candidate before the one under consideration needs to
        #       be the same
        # 2. Someone faces off against c and is defeated by c
        # 3. Then we sum over the probability of every variation other than
        #       that
        seqMap = {}
        for seq in mainSeqs:
            '''
            Make a list of all the possible alternative sequences. These are 
            sequences where whoever faces off against the candidate falls in
            that position instead of the candidate under consideration. We also
            do not need to compute probabilities of alternative sequences that
            result in the same winner being chosen, because any such sequence
            will by definition have expected utility 0 regardless of its
            probability. It won't affect the results, but for the sake of
            performance we can skip those computations
            '''
            candLoc = [i for i,x in enumerate(seq) if x == c][0]
            altStart = list(seq[:candLoc])
            lastCands = seq[candLoc:]
            altSeqs = [list(i) for i in it.permutations(lastCands) \
                               if (i[-1] != seq[-1] and i[-1] != c)]
            altSeqEnds = copy.deepcopy(altSeqs)
            for i in range(len(altSeqs)):
                altSeqs[i] = altStart + altSeqs[i]
            seqMap[ListToStr(seq)] = altSeqs

            #The question of indirect pivotality is as follows. What is the
            # probability of my vote for the candidate changing the outcome
            #  from one of the main sequences to one of the alternative
            #  sequences, and what is the utility difference of that result?
            #To find this, first get the baseline starting vote total
            initVotes = GetVotes(all_ballots, cands, [])
            #Start the probability as the probability that the base
            # sequence under consideration will happen. The dropped
            # candidates are all the candidates up to the last one
            baseProb = all_drop_probs[ListToStr(seq[:-1])]
            
            #Now for each alternative sequence, find the probability of
            # that sequence occurring and if we just gave c one more vote
            for altSeq in altSeqs:
                #Because we neglect more-than-2-way ties, the opponent in
                # the pivotal event must be the candidate dropped in
                # position candLoc instead of cand
                opponent = altSeq[candLoc]
                #First requirement is that we are raising c up above the
                # candidate that would otherwise defeat them
                altOpp = altSeq[0]
                breakSwitchProb = skpmf(0, initVotes[c], initVotes[opponent])
                makeSwitchProb = skpmf(0, initVotes[c], initVotes[opponent]-1)
                #Now this candidate gets dropped
                dropped = [altOpp]
                votes = GetVotes(all_ballots, cands, dropped)
                #Then we want to know the probability of this drop sequence
                # given the drop(s) that have already happened. This is all
                # of the candidates in this drop order except the last,
                # concatenated to the end of the preceding drop(s)
                newDropSeq = altSeq[:-1]
                newDropProb = all_drop_probs[ListToStr(newDropSeq)]
                fullAltProb = baseProb * (breakSwitchProb + \
                                          0.5*makeSwitchProb) * newDropProb

                #Record the result
                indirPivotProbs[ballotLabel] += fullAltProb
                
                if seq == switch_to_print["base"] and \
                   altSeq == switch_to_print["alt"]:
                    print(f"{c} defeats {opponent}, switching from {seq}"+\
                          f" to {altSeq} with probability {fullAltProb}")
    return(indirPivotProbs)

def CalcSmdPivotProb(ballots, held_ballots, cands):
    votes = {c:0 for c in cands}
    for b in held_ballots:
        top = b[0]
        votes[top] += 1
    all_pivot_probs = {ListToStr(b):None for b in ballots}
    for b in ballots:
        top = b[0]
        others = [c for c in cands if c != top]
        pivotProb = 0
        for alt in others:
            topProb = 1
            #What happens if this candidate is the first-place competitor
            for rem in others:
                #This requires that it has more votes than any other competitor
                if alt != rem:
                    topProb *= skcdf(-1, votes[rem], votes[alt])
            breakProb = skpmf(0, votes[top], votes[alt])
            makeProb = skpmf(0, votes[top], votes[alt]-1)
            pivotProb += topProb * (breakProb + 0.5*makeProb)
        all_pivot_probs[ListToStr(b)] = copy.deepcopy(pivotProb)
    return(all_pivot_probs)

def JitterList(xVals, epsilon):
    _jitteredList = [None] * len(xVals)
    for i in range(len(_jitteredList)):
        _jitteredList[i] = xVals[i] + np.random.uniform(-epsilon,epsilon)
    return(_jitteredList)


###############################################################################
# DEFAULT SCENARIO
###############################################################################
HELD_BALLOTS = RepList(["gross", "palin", "begich", "claus", "peltola"], 2) +\
               RepList(["claus", "gross", "palin", "begich", "peltola"], 2) +\
               RepList(["begich", "gross", "palin", "claus", "peltola"], 3) +\
               RepList(["palin", "peltola", "begich", "claus", "gross"], 6) +\
               RepList(["peltola", "begich", "claus", "gross", "palin"], 12) 

print("DEFAULT SCENARIO simulation:")
IRV(cands = CANDS,
    ballots = HELD_BALLOTS)
print("\n")


###############################################################################
# DIRECTLY PIVOTAL SCENARIO
###############################################################################
DIR_BALLOT = ["peltola", "begich", "claus", "gross", "palin"]
DIR_HELD_BALLOTS = HELD_BALLOTS + [DIR_BALLOT]

print("DIRECTLY PIVOTAL simulation:")
IRV(cands = CANDS,
    ballots = DIR_HELD_BALLOTS)
print("\n")


###############################################################################
# DIRECT PIVOTALITY CALCULATION
###############################################################################
#We only care about one drop sequence. At issue is the probability that there
# is a tie for last, so we only want to consider the probability of the first
# three drops
DIR_SEQUENCE = ["gross","claus","begich"]
suppCand = "peltola"
finalOpp = "palin"

#Direct pivotal probabilities
dropSeqProb = CalcDropProbs(sequences = [DIR_SEQUENCE], \
              ballots = HELD_BALLOTS, \
              cands = CANDS)[ListToStr(DIR_SEQUENCE)]

#The candidates which still have votes are c and its direct opponent
votes = GetVotes(ballots = HELD_BALLOTS, \
                 cands = CANDS, \
                 dropped = DIR_SEQUENCE)

#Now that we have the probability of this sequence of drops occurring,
# we need only to know the probability that a vote for the candidate
# under consideration will be pivotal at the end of this chain
breakProb = dropSeqProb * skpmf(0, votes[suppCand], votes[finalOpp])
makeProb = 0.5 * dropSeqProb * skpmf(-1, votes[suppCand], votes[finalOpp])

print(f"PROBABILITY OF DIRECTLY PIVOTAL EXAMPLE: \n" +\
      f"{breakProb + makeProb} \n\n")


###############################################################################
# DOUBLE-CHECK MANUAL CALCULATION OF DIRECT PIVOT PROBABILITY USING FUNCTION
###############################################################################
print("AUTOMATIC COMPUTATION OF THE DIRECTLY PIVOTAL PROBABILITY:")
all_drop_seqs = []
for i in range(1,len(CANDS)):
    all_drop_seqs += [list(i) for i in it.permutations(CANDS, i)]
all_drop_probs = CalcDropProbs(sequences = all_drop_seqs, \
                         ballots = HELD_BALLOTS, \
                         cands = CANDS)
CalcDirectPiv(cands = CANDS,
              ballot = ["peltola", "gross", "begich"],
              all_ballots = HELD_BALLOTS,
              all_drop_probs = all_drop_probs,
              probs_to_print = [["gross", "claus", "begich", "palin"]])
print("\n")


###############################################################################
# INDIRECTLY PIVOTAL SCENARIO
###############################################################################
INDIR_BALLOT = ["gross", "begich", "claus", "peltola", "palin"]
INDIR_HELD_BALLOTS = HELD_BALLOTS + [INDIR_BALLOT]

print("INDIRECTLY PIVOTAL simulation:")
IRV(cands = CANDS,
    ballots = INDIR_HELD_BALLOTS)
print("\n")


###############################################################################
# INDIRECT PIVOTALITY CALCULATION
###############################################################################
BASE_SEQ = ["gross", "claus", "begich", "peltola", "palin"]
ALT_SEQ = ["claus", "begich", "palin", "gross", "peltola"]
suppCand = "gross"
candLoc = 0

initVotes = GetVotes(ballots = INDIR_HELD_BALLOTS, \
                     cands = CANDS, \
                     dropped = [])

baseProb = CalcDropProbs(sequences = [BASE_SEQ], \
                         ballots = HELD_BALLOTS, \
                         cands = CANDS)[ListToStr(BASE_SEQ)]
opponent = ALT_SEQ[candLoc]

breakProb = skpmf(0, initVotes[suppCand]-1, initVotes[opponent])
makeProb = skpmf(-1, initVotes[suppCand]-1, initVotes[opponent])
dropped = [opponent]

votes = GetVotes(ballots = INDIR_HELD_BALLOTS, \
                 cands = CANDS, \
                 dropped = dropped)

altDropSeq = ALT_SEQ[:-1]
altDropProb = CalcDropProbs(sequences = [ALT_SEQ], \
                            ballots = HELD_BALLOTS, \
                            cands = CANDS)[ListToStr(ALT_SEQ)]

fullAltProb = baseProb * (breakProb + 0.5*makeProb) * altDropProb

print(f"PROBABILITY OF INDIRECTLY PIVOTAL EXAMPLE: \n {fullAltProb} \n\n")


###############################################################################
# DOUBLE-CHECK MANUAL CALCULATION OF INDIRECT PIVOT PROBABILITY USING FUNCTION
###############################################################################
print("AUTOMATIC COMPUTATION OF THE INDIRECTLY PIVOTAL PROBABILITY:")

CalcIndirectPiv(cands = CANDS,
                ballot = ["gross", "begich", "claus", "peltola", "palin"],
                all_ballots = HELD_BALLOTS,
                all_drop_probs = all_drop_probs,
                switch_to_print = {"base": ["gross", "claus", "begich",
                                            "peltola", "palin"], \
                                   "alt": ["claus", "begich", "palin", "gross",
                                           "peltola"]})


###############################################################################
# SIMULATE PIVOTALITY
###############################################################################
print("\n\n SIMULATING PIVOT PROBABILITIES \n\n")

POSSIBLE_CANDS = ["A", "B", "C", "D", "E", "F"]
cand_nums = [3, 4, 5]

dirLessSmdTotal = 0
np.random.seed(20221003)
powerTau = 0.5

run_types = ['uniform','power']
for runType in run_types:
    probs_to_plot = {n:{'irv': [], 'smd': []} for n in cand_nums}
    if runType == 'power':
        PLOT_EPSILON = 0.15
        YMIN = -15
        YMAX = 0
    elif runType == 'uniform':
        PLOT_EPSILON = 0.15
        YMIN = 0
        YMAX = 0.12

    for run in range(NUM_RUNS):
        if SEED is not None:
            np.random.seed(SEED)
        else:
            SEED = np.random.randint(2**30)
            np.random.seed(SEED)
            print(f"Seed for run {run}: {SEED}")
            SEED = None
        
        runStart = time.time()
        for candNum in cand_nums:
            CANDS = POSSIBLE_CANDS[:candNum]
            DEFAULT_BALLOT = copy.deepcopy(CANDS)
            NUM_BALLOTS = 100
            HELD_BALLOTS = [None] * NUM_BALLOTS
            if runType == 'uniform':
                for i in range(NUM_BALLOTS):
                    np.random.shuffle(DEFAULT_BALLOT)
                    HELD_BALLOTS[i] = copy.deepcopy(DEFAULT_BALLOT)
            elif runType == 'power':
                for i in range(NUM_BALLOTS):
                    if i <= int(np.floor(NUM_BALLOTS*powerTau)):
                        np.random.shuffle(DEFAULT_BALLOT)
                    HELD_BALLOTS[i] = copy.deepcopy(DEFAULT_BALLOT)
            #Get all drop probabilities
            all_drop_seqs = []
            for i in range(1,len(CANDS)):
                all_drop_seqs += [list(i) for i in it.permutations(CANDS, i)]
            
            all_drop_probs = CalcDropProbs(sequences = all_drop_seqs,\
                                           ballots = HELD_BALLOTS,\
                                           cands = CANDS)
            #print(all_drop_probs) 
            diffThisRun = 0
            irvPivotsThisRun = []
            smdPivotsThisRun = []
            for ballot in [i for i in it.permutations(DEFAULT_BALLOT)]:
                dir_pivot_probs = CalcDirectPiv(cands = CANDS,
                                            ballot = ballot,
                                            all_ballots = HELD_BALLOTS,
                                            all_drop_probs = all_drop_probs,
                                            probs_to_print = [])
                smd_pivot_probs = CalcSmdPivotProb([ballot],HELD_BALLOTS,CANDS)
                indir_pivot_probs = CalcIndirectPiv(cands = CANDS,
                                            ballot = ballot,
                                            all_ballots = HELD_BALLOTS,
                                            all_drop_probs = all_drop_probs,
                                            switch_to_print = NO_SWITCHES)
                diff = dir_pivot_probs[ListToStr(ballot)] - \
                       smd_pivot_probs[ListToStr(ballot)]
                diffThisRun += diff
                irvPivotsThisRun.append(dir_pivot_probs[ListToStr(ballot)] + \
                                        indir_pivot_probs[ListToStr(ballot)])
                smdPivotsThisRun.append(smd_pivot_probs[ListToStr(ballot)])
                if VERBOSE:
                    print(f"DIRECT PIVOT PROBABILITY: {dir_pivot_probs}")
                    print(f"INDIRECT PIVOT PROBABILITY: {indir_pivot_probs}")
                    print(f"SMD PIVOT PROBABILITY: {smd_pivot_probs}")
                    print("\n")
            dirLessSmdTotal += diffThisRun
            probs_to_plot[candNum]['irv'].append(np.mean(irvPivotsThisRun))
            probs_to_plot[candNum]['smd'].append(np.mean(smdPivotsThisRun))
            if ERROR_CHECKS:
                for uNums in range(1,len(CANDS)):
                    dropSeqSum = sum(all_drop_probs[k] for k in \
                        [_ for _ in all_drop_probs.keys() if len(_) == uNums])
                    print(\
                      f"Probability of ANY drop sequence of length {uNums}:"+\
                      f" {dropSeqSum} with {candNum} candidates.")
                print("\n")
            with open(f"probs_to_plot_{runType}.txt", "w") as f:
                f.write(str(probs_to_plot))
            
        elapsed = round((time.time() - runStart)/60,2)
        print(f"Finished run {run+1} of {NUM_RUNS} after {elapsed} minutes")
    fig, ax = plt.subplots(figsize=(9,7))
    ax.set_xlim(min(cand_nums)-0.5, max(cand_nums)+0.5)
    ax.set_ylim(YMIN, YMAX)

    if runType == 'uniform':
        ax.set_title('Uniform support')
        ax.set_ylabel("Pivotal probability", labelpad=15)
    elif runType == 'power':
        ax.set_title('Power law support')
        ax.set_ylabel("Pivotal probability (Log$_{10}$)", labelpad = 15)

    ax.set_xlabel("Number of candidates", labelpad=10)
    xTickLocs = []
    xTickLabels = []
    for candNum in cand_nums:
        xTickLocs.append(candNum)
        xTickLabels.append(str(candNum))
    plt.xticks(xTickLocs)
    plt.gca().axes.xaxis.set_ticklabels(xTickLabels)
    plt.tight_layout()
    IRVCOL = (237/255, 6/255, 119/255, 0.5)
    SMDPCOL = (0/255, 158/255, 219/255, 0.5)
    MARKERSIZE = 10
    
    plt.axvline(x=3.5)
    plt.axvline(x=4.5)

    for candNum in cand_nums:
        x = list(np.repeat(candNum, len(probs_to_plot[candNum]['smd'])))
        originalX = copy.deepcopy(x)
        
        if runType == 'uniform':
            x = [_ - 0.2 for _ in JitterList(originalX, PLOT_EPSILON)]
            ax.plot(x, probs_to_plot[candNum]['smd'], marker='s', 
                    linestyle='None', markersize=MARKERSIZE, 
                    color=SMDPCOL, markeredgewidth=1, markeredgecolor="black")
            x = [_ + 0.2 for _ in JitterList(originalX, PLOT_EPSILON)]
            ax.plot(x, probs_to_plot[candNum]['irv'], marker='o', 
                    linestyle='None', markersize=MARKERSIZE, 
                    color=IRVCOL, markeredgewidth=1, markeredgecolor="black")
            currSmdMedian = np.median(probs_to_plot[candNum]['smd'])
            currIrvMedian = np.median(probs_to_plot[candNum]['irv'])
        elif runType == 'power':
            x = JitterList(originalX, PLOT_EPSILON)
            valsToMed = {'smd': [], 'irv': []}
            for i in range(len(x)):
                nextVal = mpmath.log(probs_to_plot[candNum]['smd'][i], 10)
                valsToMed['smd'].append(nextVal)
                pair = (x[i], nextVal)
                ax.plot(pair[0]-0.2, pair[1], marker='s', linestyle='None',
                        markersize=MARKERSIZE, color=SMDPCOL,
                        markeredgewidth=1, markeredgecolor="black")
            x = JitterList(originalX, PLOT_EPSILON)
            for i in range(len(x)):
                nextVal = mpmath.log(probs_to_plot[candNum]['irv'][i], 10)
                valsToMed['irv'].append(nextVal)
                pair = (x[i], nextVal)
                ax.plot(pair[0]+0.2, pair[1], marker='o', linestyle='None', 
                        markersize=MARKERSIZE, color=IRVCOL, 
                        markeredgewidth=1, markeredgecolor="black") 
            currSmdMedian = np.median(valsToMed['smd'])
            currIrvMedian = np.median(valsToMed['irv'])
        
        plt.hlines(y=currSmdMedian, xmin=candNum-0.5, xmax=candNum, 
                   colors=SMDPCOL, linewidth=3)
        plt.hlines(y=currIrvMedian, xmin=candNum, xmax=candNum+0.5,
                   colors=IRVCOL, linewidth=3)
        
        ax.legend(handles = [Line2D([0], [0], marker = 's', 
                  markerfacecolor = SMDPCOL, markersize=MARKERSIZE, 
                  markeredgewidth=1, markeredgecolor="black", 
                  label = "SMDP", lw=0),
                  Line2D([0], [0], marker = 'o', markersize=MARKERSIZE,
                         markerfacecolor = IRVCOL, markeredgewidth=1, 
                         markeredgecolor="black", label = "IRV", lw=0)],
                  loc="upper right",
                  borderpad=0.5,
                  labelspacing=0.5)
    plt.savefig('../results/simulated_probabilities_'+runType+'.png')
